close all;
clear;
%% 参数设置
T = 1e-3;               % 总时长
fs = 26e6;              % 采样率26MHz
gps_fc = 6.5e6;         % 中频频率
gps_cacode_fc = 1023e3; % CA码码率
gps_f_err = 2.0e2;      % 生成GPS信号时时的频率偏移

%% 读取信号
DataFilePath = "./gps_data.mat";
N = T/(1/fs);   % 采样点数

gps_data = load(DataFilePath);
gps_sig = gps_data.gps_sig;

%% 频谱
[sig_f, sig_p] = fft_calc(gps_sig, fs, 26000);

figure(1);
subplot(211);plot(gps_sig);xlabel('N');ylabel('A');title('时域波形');
subplot(212);plot(sig_f, sig_p);xlabel('f (Hz)');ylabel('P1(f)');title('调制后频谱');

%% 本地伪码/载波
local_NCO_sin = zeros(1, 26000);
local_NCO_cos = zeros(1, 26000);

for i = 1:26000
    cnt = ceil(i * gps_cacode_fc/fs);   % 计算CA码下标
    local_NCO_sin(i) = sin(2*pi*gps_fc*(i*T/26000));
    local_NCO_cos(i) = cos(2*pi*gps_fc*(i*T/26000));
end

%% 混频
gps_data_i = gps_sig.*local_NCO_cos;
gps_data_q = gps_sig.*local_NCO_sin;

%% 低通滤波
filter_i = decm_filter();    % 加载滤波器
[bi, ai] = tf(filter_i);    % 传递函数
gps_data_i_lpf = filter(bi, ai, gps_data_i);


filter_q = decm_filter();    % 加载滤波器
[bq, aq] = tf(filter_q);    % 传递函数
gps_data_q_lpf = filter(bq, aq, gps_data_q);

%% 降采样
gps_data_i_decm = zeros(1, 2046);
gps_data_q_decm = zeros(1, 2046);
temp = 0;
for i = 1:26000
    cnt = ceil(i * gps_cacode_fc*2/fs);   % 计算CA码下标
    if temp == cnt-1
        gps_data_i_decm(cnt) = gps_data_i_lpf(i);
        gps_data_q_decm(cnt) = gps_data_q_lpf(i);
        temp = cnt;
    end
end


%% 混频、滤波、抽样后频谱
[gps_sig_f, gps_sig_p] = fft_calc(gps_sig, fs, 26000);
[gps_data_iq_f, gps_data_iq_p] = fft_calc((gps_data_i + gps_data_q*1i), fs, 26000);
[gps_data_iq_lpf_f, gps_data_iq_lpf_p] = fft_calc((gps_data_i_lpf + gps_data_q_lpf*1i), fs, 26000);
[gps_data_iq_decm_f, gps_data_iq_decm_p] = fft_calc((gps_data_i_decm + gps_data_q_decm*1i), 2046000, 2046);

figure(2);
subplot(411);plot(gps_sig_f, gps_sig_p);xlabel('f (Hz)');ylabel('P1(f)');title('混频前频谱');
subplot(412);plot(gps_data_iq_f, gps_data_iq_p);xlabel('f (Hz)');ylabel('P1(f)');title('混频后频谱');
subplot(413);plot(gps_data_iq_lpf_f, gps_data_iq_lpf_p);xlabel('f (Hz)');ylabel('P1(f)');title('滤波后频谱');
subplot(414);plot(gps_data_iq_decm_f, gps_data_iq_decm_p);xlabel('f (Hz)');ylabel('P1(f)');title('抽样后频谱');
figure(3);
subplot(411);plot(gps_sig(1:26000));xlabel('N');ylabel('A');title('混频前波形');
subplot(412);plot(gps_data_i(1:26000));xlabel('N');ylabel('A');title('混频后波形');
subplot(413);plot(gps_data_i_lpf(1:26000));xlabel('N');ylabel('A');title('滤波后波形');
subplot(414);plot(gps_data_i_decm);xlabel('N');ylabel('A');title('抽样后波形');

%% 相干积分
fft_N = 256;
groupSize = 22;
groupNum = 2046/groupSize;

gps_data_res = zeros(1023, 2048);
gps_cacode_2046K = zeros(1, 2046);
temp = GPS_L1_CACode([2 3]);
gps_data_sum = zeros(1023, groupNum);

for x = 1:1023
    for i = 1:1023
        gps_cacode_2046K(i*2-1: i*2) = [temp(i) temp(i)];
    end

    for i = 1:groupNum
        temp_i = gps_data_i_decm((i-1)*groupSize+1:i*groupSize).*gps_cacode_2046K((i-1)*groupSize+1:i*groupSize);
        temp_q = gps_data_q_decm((i-1)*groupSize+1:i*groupSize).*gps_cacode_2046K((i-1)*groupSize+1:i*groupSize);
        gps_data_sum(x, i) = sum(temp_i + temp_q*1i);
    end
    temp_ = fft(gps_data_sum(x, :), fft_N);
    gps_data_res(x, 1:fft_N) = abs(temp_);

    temp(1:1023) = [temp(2:1023) temp(1)];

    disp(x);
end

x_coord = 2046e3/groupSize * (1 - (1:fft_N)./fft_N);
y_coord = 1:1023;
figure(4);mesh(x_coord, y_coord, gps_data_res(:, 1:fft_N));ylabel('PH');xlabel('f(Hz)');

    
% x --> j
% y --> i
max_x = 0;
max_y = 0;
max_A = 0;

for i = 1:1023
    for j = 1:fft_N
        if gps_data_res(i, j) > max_A
            max_x = j;
            max_y = i;
            max_A = gps_data_res(i, j);
%             disp(i);
        end
    end
end

disp("(频率)最高点X");
disp(max_x);
disp("对应频率：");
disp(2046e3/groupSize * (1 - max_x/fft_N));
disp("(相位)最高点Y");
disp(max_y);







































